Unsupervised Classification Methods for lidar-based tree crown deliniation

tile=readLAS("../tests/data/NEON_D03_OSBS_DP1_405000_3276000_classified_point_cloud.laz")

tile<-tile %>% lasfilter(Z < 100)
tile<-tile %>% lasfilter(Z > 0)

plot(tile, color="Z", colorPalette = heat.colors(50), trim=0.95)

#remove 
summary(tile)
#table(tile@data$Classification)
#ggplot(tile@data,aes(x=as.factor(Classification),y=Z)) + geom_boxplot()

#add a point index
tile@data$PointID<-1:nrow(tile@data)

#Lastrees updates the cloud inplace, save data in seperate object
ptlist<-list()

#ground model
ground_model(las)

plot(li2012,color="treeID",colorPalette = sample(pastel.colors(500)), size = 1)

#canopy model
chm=canopy_model(las)

#Toy file for testing
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
tile = readLAS(LASfile, select = "xyz", filter = "-drop_z_below 0")
col = pastel.colors(200)

#add a point index
tile@data$PointID<-1:nrow(tile@data)

#Lastrees updates the cloud inplace, save data in seperate object
ptlist<-list()

ground_model(tile)

chm=canopy_model(tile)

plot(tile,color="Classification")

You must enable Javascript to view this page properly.

Watershed

system.time(watershed<-segment_ITC(tile,algorithm = "watershed",chm=chm))

##    user  system elapsed 
##   5.438   0.098  19.046

ptlist[["watershed"]]<-watershed %>% lasfilter(!is.na(treeID)) %>% .@data
plot(watershed,color="treeID",colorPalette = sample(pastel.colors(500)), size = 1)

Li 2012

system.time(li2012<-segment_ITC(tile,algorithm = "li2012",chm=chm))

##    user  system elapsed 
##   0.284   0.025   1.217

ptlist[["li2012"]]<-li2012 %>% lasfilter(!is.na(treeID)) %>% .@data
plot(li2012,color="treeID",colorPalette = sample(pastel.colors(500)), size = 1)

Dalponte 2016

system.time(dalponte2016<-segment_ITC(tile,algorithm = "dalponte2016",chm=chm))

##    user  system elapsed 
##   4.301   0.051  17.771

ptlist[["dalponte2016"]]<-dalponte2016 %>% lasfilter(!is.na(treeID)) %>% .@data

Silva 2016

system.time(silva2016<-segment_ITC(tile,algorithm = "silva2016",chm=chm))

##    user  system elapsed 
##   4.268   0.054  17.102

ptlist[["silva2016"]]<-silva2016 %>% lasfilter(!is.na(treeID)) %>% .@data 

Comparison

How many tree objects

library(pander)
df<-data.frame(Algorthm=c("watershed","li2012","Dalponte2016","Silva2016"),Total_Trees=c( max(ptlist[["watershed"]]$treeID,na.rm=T),max(ptlist[["li2012"]]$treeID,na.rm=T), max(ptlist[["dalponte2016"]]$treeID,na.rm=T), max(ptlist[["silva2016"]]$treeID,na.rm=T)))
pandoc.table(df,style="rmarkdown")

## 
## 
## |   Algorthm   | Total_Trees |
## |:------------:|:-----------:|
## |  watershed   |     76      |
## |    li2012    |     378     |
## | Dalponte2016 |     111     |
## |  Silva2016   |     111     |

Consensus

pdf<-melt(ptlist,id.vars=colnames(ptlist[["silva2016"]]))
head(pdf)

##          X       Y      Z PointID Classification  Zref treeID        L1
## 1 481352.0 3813013  0.000      23              2  0.33      2 watershed
## 2 481352.0 3813013  0.000      24              2  0.60      2 watershed
## 3 481352.3 3813015 14.289      25              0 15.18      2 watershed
## 4 481352.3 3813015 14.835      26              0 15.49      2 watershed
## 5 481352.2 3813014 14.595      27              0 15.02      2 watershed
## 6 481352.2 3813015 14.409      28              0 14.97      2 watershed

Ground versus Trees

Since the input of segmentation methods differ among algorithms, there may be different thresholds for what constitutes a tree. Let’s remove points that only appear in less than 3 methods.

#How many methods does each point appear in.
table(table(pdf$PointID))

## 
##     1     2     3     4 
##  5342  2070  2976 35239

points_to_remove<-pdf %>% group_by(PointID) %>% summarize(n=n()) %>% arrange(n) %>% filter(n < 3) %>% .$PointID

pdf<-pdf %>% filter(!PointID %in% points_to_remove)

Did that improve the similarity of # of tree objects?

pdf %>% group_by(L1) %>% summarize(Total_Trees=length(unique(treeID)))

## # A tibble: 4 x 2
##   L1           Total_Trees
##   <chr>              <int>
## 1 dalponte2016         111
## 2 li2012               356
## 3 silva2016            111
## 4 watershed             76

Majority Rule

#Create point ID lookup table
res<-dcast(pdf,PointID~L1,value.var = "treeID")

idframe<-res %>% add_rownames() %>% select(rowname,PointID)
head(res<-res %>% select(-PointID))

##   dalponte2016 li2012 silva2016 watershed
## 1          110     35         5         2
## 2          110     35        17         2
## 3          111     35         5         2
## 4          111     35         5         2
## 5          110     35         5         2
## 6          111     35         5         2

system.time(consensus<-majority_voting(res, is.relabelled = FALSE))

##    user  system elapsed 
##   6.648   0.105  20.573

#reassign to pointID
consensus_frame<-data.frame(rowname=rownames(res),consensus=consensus) %>% inner_join(idframe) %>% select(PointID,consensus)

#merge with the original tile
tile@data<-merge(tile@data,consensus_frame,all=T)
plot(tile,color="consensus",colorPalette = sample(pastel.colors(500)), size = 1)

You must enable Javascript to view this page properly.

Next steps

• Tree Extractor (probably in python)

but let’s just look at some individual trees. The question is, can we classify them based on the lidar cloud? Do we need to drape on the RGB data.

ind_trees= split(tile@data, tile@data$treeID)
ind_trees = lapply(ind_trees, LAS, header = tile@header)

plot(ind_trees[[20]],color="treeID",bg="grey90")

You must enable Javascript to view this page properly.

plot(ind_trees[[30]],color="treeID",bg="grey90")

You must enable Javascript to view this page properly.

plot(ind_trees[[40]],color="treeID",bg="grey90")

You must enable Javascript to view this page properly.

Clearly, more work to be done. Especially in reference to the ground model.

#save.image("ITCs.RData")